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Abstract. In this paper, we construct new finite element methods for the approximation 
of the equations of linear elasticity in three space dimensions that produce direct approxima- 
tions to both stresses and displacements. The methods are based on a modified form of the 
Hcllinger-Reissner variational principle that only weakly imposes the symmetry condition 
on the stresses. Although this approach has been previously used by a number of authors, 
a key new ingredient here is a constructive derivation of the elasticity complex starting 
from the de Rham complex. By mimicking this construction in the discrete case, we derive 
new mixed finite elements for elasticity in a systematic manner from known discretizations 
of the de Rham complex. These elements appear to be simpler than the ones previously 
derived. For example, we construct stable discretizations which use only piecewise linear 
elements to approximate the stress field and piecewise constant functions to approximate 
the displacement field. 



1. Introduction 

The equations of linear elasticity can be written as a system of equations of the form 

(1.1) Aa = eu, div a = f in Q. 

Here the unknowns a and u denote the stress and displacement fields engendered by a body 
force / acting on a linearly elastic body which occupies a region (let 3 . Then a takes values 
in the space § := of symmetric matrices and u takes values in V := R 3 . The differential 
operator e is the symmetric part of the gradient, the div operator is applied row- wise to a 
matrix, and the compliance tensor A = A(x) : 8 — > § is a bounded and symmetric, uniformly 
positive definite operator reflecting the properties of the body. If the body is clamped on the 
boundary dQ of Q, then the proper boundary condition for the system (II. ip is u = on dQ. 
For simplicity, this boundary condition will be assumed here. The issues that arise when 
other boundary conditions are assumed (e.g., the case of pure traction boundary conditions 
an = g) are discussed in [9]. 

The pair (a, u) can alternatively be characterized as the unique critical point of the 
Hellinger-Reissner functional 

(1.2) J(r,v) = / (-At : r + div r • v — f • v) dx. 
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The critical point is sought among all r e H(div, Q; §), the space of square-integrable sym- 
metric matrix fields with square-integrable divergence, and all v € L 2 (Q;Y), the space of 
square-integrable vector fields. Equivalently, ((T,u) G if(div, fi;S) x L 2 (Q;Y) is the unique 
solution to the following weak formulation of the system (11. ip : 



A mixed finite element method determines an approximate stress field ah and an approx- 
imate displacement field Uh as the critical point of J over E^ x Vh where E fe C if (div, Q; S) 
and Vh C L 2 (Q;Y) are suitable piecewise polynomial subspaces. Equivalently, the pair 
(<7h, u h) G E ft x Vh is determined by the weak formulation (ll.3p . with the test space re- 
stricted to E/, x V),. As is well known, the subspaces E^ and Vh cannot be chosen arbitrarily. 
To ensure that a unique critical point exists and that it provides a good approximation of 
the true solution, they must satisfy the stability conditions from Brezzi's theory of mixed 
methods [HICES]. 

Despite four decades of effort, no stable simple mixed finite element spaces for elasticity 
have been constructed. For the corresponding problem in two space dimensions, stable finite 
elements were presented in [10]. For the lowest order element, the space E^ is composed 
of piecewise cubic functions, with 24 degrees of freedom per triangle, while the space Vh 
consists of piecewise linear functions. Another approach which has been discussed in the 
two dimensional case is the use of composite elements, in which Vh consists of piecewise 
polynomials with respect to one triangulation of the domain, while E ft consists of piecewise 
polynomials with respect to a different, more refined, triangulation [5J, E31 EI] • In three 
dimensions, a partial analogue of the element in [10] has been proposed and shown to be 
stable in [I]. This element uses piecewise quartic stresses with 162 degrees of freedom per 
tetrahedron, and piecewise linear displacements. 

Because of the lack of suitable mixed elasticity elements, several authors have resorted to 
the use of Lagrangian functionals which are modifications of the Hellinger-Reissner functional 
given above [21 HI El [271 12S I2S1 EQ] , in which the symmetry of the stress tensor is enforced 
only weakly or abandoned altogether. In order to discuss these methods, we consider the 
compliance tensor A(x) as a symmetric and positive definite operator mapping M into M, 
where M is the space of 3 x 3 matrices. In the isotropic case, for example, the mapping 
a 1— > Aa has the form 



where X(x),/j,(x) are positive scalar coefficients, the Lame coefficients. A modification of 
the variational principle discussed above is obtained if we consider the extended Hellinger- 
Reissner functional 






J n div a ■ vdx 




(1.4) 




over the space H(div, 0;M) x L 2 (Q;Y) x L 2 (Q;K), where K denotes the space of skew 
symmetric matrices. We note that the symmetry condition for the space of matrix fields is 
now enforced through the introduction of the Lagrange multiplier, q. A critical point (a, u,p) 
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of the functional J e is characterized as the unique solution of the system 



(1.5) 



f Q (Aa : r + div t ■ u + r : p) dx 
J n div a ■ vdx 
J Q a: qdx 



0, r G H(div, Q; M) 

J n f-vdx, veL 2 (n;Y), 
0, q G L 2 {Vl; K). 



Clearly, if (a, u,p) is a solution of this system, then a is symmetric, i.e., a G //(div, f2; §), and 
therefore the pair (a, it) G if (div, f2; S) x L 2 (f2; V) solves the corresponding system (jl.3p . On 
the other hand, if (it, p) solves (ll.3p . then u G i/ 1 (f2; V) and, if we set p to the skew-symmetric 
part of gradw, then (cr,u,p) solves (11.51) . In this respect, the two systems (II. 3p and (11.51) are 
equivalent. However, the extended system (I1.5P leads to new possibilities for discretization. 
Assume that we choose finite element spaces Ef, x 14 x Qh C i?(div, f2;M) x L 2 (Q;Y) x 
L 2 (Q; IK) and consider a discrete system corresponding to (11.51) . If (<T/j, u h ,p h ) G T, h x V h x Q h 
is a discrete solution, then ah will not necessary inherit the symmetry property of a. Instead, 
ah will satisfy the weak symmetry condition 



Therefore, these solutions will in general not correspond to solutions of the discrete system 
obtained from (11.31) . 

Discretizations based on the system (II. 5p will be referred to as mixed finite element meth- 
ods with weakly imposed symmetry. Such discretizations were already introduced by Fraejis 
de Veubeke in [21] and further developed in [1]. In particular, the so-called PEERS element 
proposed in [3] for the corresponding problem in two space dimensions used a combination of 
piecewise linear functions and cubic bubble functions, with respect to a triangulation of the 
domain, to approximate the stress a, piecewise constants to approximate the displacements, 
and continuous piecewise linear functions to approximate the Lagrange multiplier p. Prior 
to the PEERS paper, Amara and Thomas [2] developed methods with weakly imposed sym- 
metry using a dual hybrid approach. The lowest order method they discussed approximates 
the stresses with quadratic polynomials plus bubble functions and the multiplier by discon- 
tinuous constant or linear polynomials. The displacements are approximated on boundary 
edges by linear functions. Generalizations of the idea of weakly imposed symmetry to other 
triangular elements, rectangular elements, and three space dimensions were developed in 
[28] . [29] . [30] and [21]. In [29], a family of elements is developed in both two and three 
dimensions. The lowest order element in the family uses quadratics plus the curls of quar- 
tic bubble functions in two dimensions or quintic bubble functions in three dimensions to 
approximate the stresses, discontinuous linears to approximate the displacements, and dis- 
continuous quadratics to approximate the multiplier. In addition, a lower order method is 
introduced that approximates the stress by piecewise linear functions augmented by the curls 
of cubic bubble functions plus a cubic bubble times the gradient of local rigid motions. The 
multiplier is approximated by discontinuous piecewise linear functions and the displacement 
by local rigid motions. Morley [23] extends PEERS to a family of triangular elements, to 
rectangular elements, and to three dimensions. In addition, the multiplier is approximated 
by nonconforming rather than continuous piecewise polynomials. 




for all q G Qh- 
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There is a close connection between mixed finite elements for linear elasticity and dis- 
cretization of an associated differential complex, the elasticity complex, which will be in- 
troduced in § 3 below. In fact, the importance of this complex was already recognized in 
[TO] , where mixed methods for elasticity in two space dimensions were discussed. The new 
ingredient here is that we utilize a constructive derivation of the elasticity complex starting 
from the de Rham complex. This construction is described in Eastwood [18] and is based 
on the the Bernstein-Gelfand-Gelfand resolution, cf. [UJ and also [14]. By mimicking the 
construction in the discrete case, we are able to derive new mixed finite elements for elas- 
ticity in a systematic manner from known discretizations of the de Rham complex. As a 
result, we can construct new elements both in two and three space dimensions which are 
significantly simpler than those derived previously. For example, we will construct stable 
discretizations of the system (jl.5p which only use piecewise linear and piecewise constant 
functions, as illustrated in the figure below. For simplicity, the entire discussion of the 
present paper will be given in the three dimensional case. A detailed discussion in two space 
dimensions can be found in [Sj. Besides the methods discussed here, we note that by slightly 
generalizing the approach of this paper, one can also analyze some of the previously known 
methods mentioned above that are also based on the weak symmetry formulation (see [Jj5] 
for details). 



Figure 1. Elements for the stress, displacement, and multiplier in the lowest 
order case in two dimensions and three dimensions. 

An alternative approach to construct finite element methods for linear elasticity is to 
consider a pure displacement formulation. Since the coefficient A in (11.11) is invertible, the 
stress a can be eliminated using the first equation in (11.11) . the stress-strain relation. This 
leads to the second order equation 



for the displacement u. A weak solution of this equation can be characterized as the global 
minimizer of the energy functional 



over the Sobolev space Hq(Q; V). Here Hq(Q;Y) denotes the space of all square integrable 
vector fields on Q, with square integrable derivatives, and which vanish on the boundary dQ. 
A finite element approach based on this formulation, where we seek a minimum over a finite 
element subspace of Hq(Q;Y) is standard and discussed in textbooks, (e.g., PS]). However, 
for more general models, arising, for example, in viscoelasticity and plasticity (cf. |15j). 
the stress-strain relation is not local and an elimination of the stress a is impossible. For 
such models, a pure displacement model is excluded, and a mixed approach seems to be 





div A eu = f in f2 
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an obvious alternative. The construction of stable mixed elements for linear elasticity is an 
important step in the construction of mixed methods for these more complicated models. 
Another advantage of the mixed approach is that we automatically obtain schemes which 
are uniformly stable in the incompressible limit, i.e., as the Lame parameter A tends to 
infinity. Since this behavior of mixed methods is well known, we will not focus further on 
this property here. A more detailed discussion in this direction can, for example, be found 
in [5]. 

An outline of the paper is as follows. In § 2, we describe the notation to be used, state 
our main result, and provide some preliminary discussion on the relation between stability 
of mixed finite element methods and discrete exact complexes. In § 3, we present two 
complexes related to the two mixed formulations of elasticity given by (11.31) and (11.51) . In 
§ 4, we introduce the framework of differential forms and show how the elasticity complex 
can be derived from the de Rham complex. In § 5, we derive discrete analogues of the 
elasticity complex beginning from discrete analogues of the de Rham complex and identify 
the required properties of the discrete spaces necessary for this construction. This procedure 
is our basic design principle. In § 6, we apply the construction of the preceding section to 
specific discrete analogues of the de Rham complex to obtain a family of discrete elasticity 
complexes. In § 7 we use this family to construct stable finite element schemes for the 
approximation of the mixed formulation of the equations of elasticity with weakly imposed 
symmetry. Finally, in § 8, we show how a slightly more complicated procedure leads to a 
simplified elasticity element. 

2. Notation, statement of main results, and preliminaries 

We begin with some basic notation and hypotheses. We continue to denote by V = M 3 
the space of 3-vectors, by M the space of 3 x 3 real matrices, and by § and IK the subspaces 
of symmetric and skew symmetric matrices, respectively. The operators sym : M — > § and 
skw : M — ► K denote the symmetric and skew symmetric parts, respectively. Note that 
an element of the space K can be identified with its axial vector in V given by the map 
vec : IK — > V: 

/ -v 3 v 2 \ I vA 
vec v 3 -v x J = \v 2 I , 
\-v 2 v x / \v 3 J 

i.e., vec _1 (f)-u; = v x w for any vectors v and w. 

We assume that Q is a domain in R 3 with boundary dQ. We shall use the standard function 
spaces, like the Lebesgue space L 2 (fl) and the Sobolev space H S (Q). For vector- valued 
functions, we include the range space in the notation following a semicolon, so L 2 (f2;X) 
denotes the space of square integrable functions mapping Q into a normed vector space X. 
The space if (div, f2; V) denotes the subspace of (vector- valued) functions in L 2 (Q; V) whose 
divergence belongs to L 2 (Q). Similarly, if (div, fl; M) denotes the subspace of (matrix- valued) 
functions in L 2 (Q; M) whose divergence (by rows) belongs to L 2 (Q; V). 

Assuming that X is an inner product space, then L 2 (Q; X) has a natural norm and inner 
product, which will be denoted by || • || and (•, •), respectively. For a Sobolev space 
ii s (fi;X), we denote the norm by || • || s and for if(div,fi;X), the norm is denoted by 



6 



DOUGLAS N. ARNOLD, RICHARD S. FALK, AND RAGNAR WINTHER 



IMIdiv := (IM| 2 + II div v || 2 ) 1//2 - The space Pfc(fi) denotes the space of polynomial functions 
on fl of total degree < k. Usually we abbreviate this to just TV 

In this paper we shall consider mixed finite element approximations derived from (II. 5p . 
These schemes take the form: 

Find (cr h ,u h ,p h ) eE h xV h X Qh such that 

f n (Aa h : r + divr • u h + r : p h ) dx = 0, r E T, h , 

(2.1) f n div a h ■ v dx = j n f-vdx, veV h , 

J n a h :qdx = 0, q G Q h , 

where now E h C if (div, O; M),V h C L 2 (fi; V) and Q h G f 2 (ft; K). 

Following the general theory of mixed finite element methods, cf. [T2J US], the stability of 
the saddle-point system (12.11) is ensured by the following conditions: 

(Al) ||t|| j iv < ci(At, r) whenever r G satisfies (divr, v) = Vi> G T4, and (r, g) = 
VgGQ fe , 

(A2) for all nonzero (v,q) G Vh X Qh, there exists nonzero r G with (divr, v) + 
(r,g)>c 2 ||r|| div (|M| + ||g||), 

where ci and c 2 are positive constants independent of h. 

The main result of this paper, given in Theorem 17. 11 is to construct a new family of stable 
finite element spaces Vh, Qh that satisfy the stability conditions (Al) and (A2). We shall 
show that for r > 0, the choices of the Nedelec second family of H(div) elements of degree 
r + 1 for Eh, cf. [26], and of discontinuous piecewise polynomials of degree r for Vh and Qh 
provide a stable finite element approximation. In contrast to the previous work described in 
the introduction, no stabilizing bubble functions are needed; nor is interelement continuity 
imposed on the multiplier. In § 8 we also discuss a somewhat simpler lowest order element 
(r = 0) in which the local stress space is a strict subspace of the full space of linear matrix 
fields. 

Our approach to the construction of stable mixed elements for elasticity is motivated by 
the success in developing stable mixed elements for steady heat conduction (i.e., the Poisson 
problem) based on discretizations of the de Rham complex. We recall (see, e.g., [7]) that 
there is a close connection between the construction of such elements and discretizations of 
the de Rham complex 

(2.2) R^C°°(ft) ^ C°°(fi;V) ^ L7°°(fi;V) ^ C°°(fi) -> 0. 

More specifically, a key to the construction and analysis of stable mixed elements is a com- 
muting diagram of the form 



C°°(fi) ^ C°°(fi;V) ^ L7°°(fi;V) ^ C°°(fi) -> 



(2.3) 



h 



W h m> U h ^ V h ^ Qh -0. 



Here, the spaces Vh C if (div) and Qh C L 2 are the finite element spaces used to discretize 
the flux and temperature fields, respectively. The spaces Uh C if (curl) and Wh C if 1 
are additional finite element spaces, which can be found for all well-known stable element 
choices. The bottom row of the diagram is a discrete de Rham complex, which is exact when 
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the de Rham complex is (i.e., when the domain is contractible). The vertical operators are 
projections determined by the natural degrees of freedom of the finite element spaces. As 
pointed out in [7], there are many such discretizations of the de Rham complex. 

A diagram analogous to (12. 3p . but with the de Rham complex replaced by the elasticity 
complex defined just below, will be crucial to our construction of stable mixed elements 
for elasticity. Discretization of the elasticity complex also gives insight into the difficulties 
of constructing finite element approximations of the mixed formulation of elasticity with 
strongly imposed symmetry, cf. [8]. 

3. The elasticity complex 

We now proceed to a description of two elasticity complexes, corresponding to strongly or 
weakly imposed symmetry of the stress tensor. In the case of strongly imposed symmetry, 
relevant to the mixed elasticity system (jl.3p . the characterization of the divergence- free 
symmetric matrix fields will be needed. In order to give such a characterization, define 
curl : C°°(f2; M) — ► C°°(Q; M) to be the differential operator defined by taking curl of each 
row of the matrix. Then define a second order differential operator J : C°°(Q; S) — ► C°°(f2; S) 
by 

(3.1) Jr = curl(curlr) T , rGC°°(fi;§). 

It is easy to check that div o J = and that J o e = 0. In other words, 

(3.2) T^C°°{Y) A C°°(§) 4 C°°(§) C°°(V) -> 0, 

is a complex. Here the dependence of the domain Q is suppressed, i.e., C°°(§) = C°°(Q;S), 
and T = T(fi) denotes the six dimensional space of infinitesimal rigid motions on Q, i.e., 
functions of the form x \— > a + Bx with a G V and Bel. In fact, when Q is contractible, 
then (13.21) is an exact sequence, a fact which will follow from the discussion below. The 
complex (13. 2p will be referred to as the elasticity complex. 

A natural approach to the construction of stable mixed finite elements for elasticity would 
be to extend the complex (13.21) to a complete commuting diagram of the form (12. 3p . where 
(13. 2p is the top row and the bottom row is a discrete analogue. However, due to the pointwise 
symmetry requirement on the discrete stresses, this construction requires piecewise polyno- 
mials of high order. For the corresponding problem in two space dimensions, such a complex 
was proposed in [10] with a piecewise cubic stress space, cf. also [5J. An analogous complex 
was derived in the three dimensional case in [5]. It uses a piecewise quartic space, with 162 
degrees of freedom on each tetrahedron for the stresses. 

We consider the formulation based on weakly imposed symmetry of the stress tensor, i.e., 
the mixed system (11.51) . Then the relevant complex is, instead of (13.21) . 

(3.3) T'^r(VxI) (grad '~ 7) > C°°(M) 4 C°°(M) (div,skw)T > r(VxK) -> 0. 
Here, 

T' = { (v, gradi>) | v G T }, 

and J : C°°(fi; M) -> C°°(fi; M) denotes the extension of the operator defined on C°°(fi; S) 
by (13. ip such that J = on C°°(f2; K). We remark that J may be written 

(3.4) Jr = curl S -1 curl r, 



cS 
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where E : 
(3.5) 



is the algebraic operator 



V = V 



tr(/i)<y, 



with 5 the identity matrix. Indeed, if r is symmetric, then curl r is trace free, and therefore 
the definition H3.4[) reduces to (13.11) on C°°(Q; §). On the other hand, if r is skew with axial 
vector u, then curlr = — Sgradti, and so curlS -1 curlr = 0. 

Observe that there is a close connection between (13.21) and (13. 31) . In fact, (13.21) can be 
derived from (13.31) by performing a projection step. To see this, consider the diagram 



r^c°°(v x k) 



(grad,-I) 



C°°(M) 



(3-6) 



7T1 



7T 2 



(diVjSkw) 2 



tm c°°(v) a c°°(S) a c°°(§) iA c°°(v) 

where the projection operators are defined by 

n (u,q) = u, 7Ti(cr) = 7r 2 (cr) = sym(a), tt 3 (u, q) = u - div q. 

We may identify C°°(V) with a subspace of C 00 (V x K), namely, 

{(u,q) : u G C°°(¥),q = skw(gradw)}. 

Under this identification, T C l7°°(V) corresponds to V C C°°(V x I 
C°° (V) on the right with a different subspace of C°° (VxK), namely, 

{(u,q): ueC°°(Y),q = 0}. 

With these identifications, the bottom row is a subcomplex of the top row, and the operators 
iTk are all projections. Furthermore, the diagram commutes. It follows easily that the 
exactness of the upper row implies exactness of the bottom row. 

In the next section, we shall discuss these complexes further. In particular, we show the 
elasticity complex with weakly imposed symmetry, i.e., (13.31) . follows from the de Rham 
complex (12. 2p . Hence, as a consequence of the discussion above, both (13. 2p and (13.31) will 
follow from (JO}. 



C7°°(V x 







0. 



We identify the 



4. From the de Rham to the elasticity complex 

The main purpose of this section is to demonstrate the connection between the de Rham 
complex (12.21) and the elasticity complex (13.31) . In particular, we show that whenever (12. 2p is 
exact, (13. 3p is exact. This section serves as an introduction to a corresponding construction 
of a discrete elasticity complex, to be given in the next section. In the following section, the 
discrete complex will be used to construct stable finite elements for the system (11.51) . 

The de Rham complex (12. 2p is most clearly stated in terms of differential forms. Here we 
briefly recall the definitions and properties we will need. We use a completely coordinate-free 
approach. For a slightly more expanded discussion and the expressions in coordinates see, 
e.g., § 4]. We let A fc denote the space of smooth differential fc-forms on f2, i.e. A = 
A k (tt) = C°°(tt; Alt fc V), where Alt fe V denotes the vector space of alternating fc-linear maps 
on V. If u) G A k we let u> x G Alt fc V denote u evaluated at x, i.e., we use subscripts to indicate 
the spatial dependence. 
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Using the inner product on Alt fe V inherited from the inner product on V (see equation 
(4.1) of [3 § 4]), we may also define the Hilbert space L 2 A fc (fi) = L 2 (fi;Alt fc V) of square 
integrable differential forms with norm denoted by || • ||, and also the rath order Sobolev 
space H m A k (£l) = H m {Vt] Alt fc V), consisting of square integrable £;-forms for which the norm 

Hi m :=(x; ii^f) i/2 

\a\<m 

is finite (where the sum is over multi-indices of degree at most ra). 

Thus, 0-forms are scalar functions and 1-forms are covector fields. We will not emphasize 
the distinction between vectors and covectors, since, given the inner product in V, we may 
identify a 1-form uj with the vector field v for which u(p) = v -p, p G V. In the 3-dimensional 
case, we can identify a 2-form uj with a vector field v and a 3-form \i with a scalar field c by 

u}(p, q) — v ■ p x q, /jL(p,q,r) — c{p x q ■ r), p,q,r£Y. 

The exterior derivative d = d k : A k — > A fc+1 is defined by 

fc+i 

(4.1) du x {v u . . . , v k+ i) = ^{-l) 3+1 d Vj u x {v u v k+1 ), to G A k , v 1 , . . . , v k+1 G V, 

3=1 

where the hat is used to indicate a suppressed argument and d v denotes the directional 
derivative in the direction of the vector v. It is useful to define 

HA k = {ueL 2 {Q; Alt fc V) | du G L 2 {Q; Alt fe+1 V) }, 

with norm given by |M|#a = IMI 2 + ll^ll 2 - Using the identifications given above, the d k 
correspond to grad, curl, and div for k = 0,1,2, respectively, and the HA k correspond to 
H\ F(curl), iy(div), and, for k = 3, L 2 . 

The de Rham complex (I2.2p can then be written 

( 4 -2) R^A° ^ A 1 ± A 2 ^ A 3 -> 0. 

It is a complex since d o d = 0. 

A differential /c-form uj on Q may be restricted to a differential /c-form on any submanifold 
M C Cl: at each point of M the restriction of u is an alternating linear form on tangent 
vectors. Moreover if dim M = k, the integral J M uj is defined. 

If X is a vector space, then A fc (X) = A fc (f2;X) refers to the fc-forms with values in X, 
i.e., A fc (X) = C°°(tt] Alt fc (V; X)), where Alt fc (V; X) are alternating A;-linear forms on V with 
values in X. Given an inner product on X, we obtain an inner product on A fc (X). Obviously 
the corresponding complex 

(4.3) X^A°(X) S A X (X) S A 2 (X) S A 3 (X) -> 0, 

is exact whenever the de Rham complex is. 

We now construct the elasticity complex as a sub complex of a complex isomorphic to the 
de Rham complex with values in the six-dimensional vector space W := K x V. First, for 
any x G M 3 we define K x : V — > K by K x v = 2skw(xf T ). We then define an operator 
K : A k {Vt-N) -> A k (tt;K) by 

(4.4) (Kw) x (vi, ...,v k ) = K x [u x (vi, . . .,v k )]. 
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Next, we define an isomorphism $ : A k (W) — > A fc (W) by 
with inverse given by 

= (u-Kn,n). 

Next, define the operator .A : A fc (W) — > A fc+1 (W) by A = $e£$ _1 . Inserting the isomorphisms 
$ in the W-valued de Rham sequence, we obtain a complex 

(4.5) $(W)^A°(W) ^ A*(W) ^ A 2 (W) ^> A 3 (W) -> 0, 

which is exact whenever the de Rham complex is. 

The operator A has a simple form. Using the definition of <&, we obtain for (u, fi) G A fe (W), 

A(u, ji) = $ o d(u; — K/i, /i) = $(<iu; — dK/i, dfj,) = (du — Sfi, dfi), 

where S = S k : A fc (V) -> A fc+1 (K), k = 0, 1,2, is given by 5 = dK-Kd. Using the definition 
(14. ip of the exterior derivative, the definition (14. 4p of K, and the Leibniz rule 

(4.6) d(u A //) = du A yu + (-l) fc ^ A d/i, a; G A fe , ^ G A £ , 
we obtain 

fc+i 

(4.7) {Su)(v x , v k+l ) = Y,{-V j+1 K Vj [u{v x , ...,v j ,... v k+l )}, u G A fc (fi; V). 

3=1 

Note that the operator S is purely algebraic, and independent of x. 
Since d 2 = we have 

= d 2 K - dKd = -{dK - Kd)d 

or 

(4.8) dS = -Sd. 
Noting that 

(S 1 fJ l )(v 1 ,v 2 ) = K Vl [fx(v 2 )] - K V2 [n(vi)} = 2skw[v 1 fi(v 2 ) T - v 2 ^{v 1 ) T ], 

fi G A^fijV), Vi,v 2 G V, 

we find, using the identity 

(4.9) a x b = —2 vec skw a6 T , 
that Si is invertible with 

(Sf^Xvi) x v 2 • w 3 = i[vec(^(v 2 ,W3)) • - vec(w(ui,U2)) • « 3 + vec(^(vi, v 3 )) ■ v 2 ], 

u G A 2 (fi; IK), Wi,w 2 ,f3 G V. 

We now define the desired subcomplex. Define 

T 1 = { (u, fi) G A x (fi; W) | du = }, T 2 = { (u, fi) G A 2 (fi; W) | = }, 
with projections tt 1 : A 1 ^; W) -> L 1 and ^ 2 : A 2 (fi; W) -> L 2 given by 
n (u, /i) = (cj, Si du), tc 2 (u, fi) = (0, /i + dS{ l u). 
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Using (14.81) . it is straightforward to check that A maps A°(W) into T 1 and T 1 into T 2 , and 
that the diagram 

$(W)^A°(W) 4. A*(W) 4 A 2 (W) 4 A 3 (W) -> 



id 


7T 1 


7T 2 


id 











$(w)^a°(w) A r 1 A r 2 A a 3 (w) -> o 

commutes, and therefore the subcomplex in the bottom row is exact when the de Rham 
complex is. This subcomplex is, essentially, the elasticity complex. Indeed, by identifying 
elements (u,fj) G T 1 with uo G A^K), and elements (Q,fi) G T 2 with /i G A 2 (V), the 
subcomplex becomes 

(4.11) $( W )^A°(KxV) ^± AHK) dl ° s " odl ; A 2 (V) A 3 (K x V) - 0. 

This complex may be identified with (13.31) . As an initial step of this identification we observe 
that the algebraic operator H : C°°(M) — > C°°(M) appearing in (13.31) via (13.41) and the 
operator Si : A : (V) — ■> A 2 (IK) are connected by the identity 

(4.12) 5 = T^SiTi, 

where Ti : C°°(M) -> A X (V) and T 2 : l7°°(M) -> A 2 (K) are given by TiF(v) = Fv and 
T 2 F(vi,v 2 ) = vec -1 F(v\ x u 2 ) for F G C°°(M). In fact, using (|4.9p . we have for any 
ui,u 2 G V 

5iT 1j P(v 1 , w 2 ) = 2skw[-Ui(Fv 2 ) T - f 2 (-Fwi) T ] 
= vec~ (i> 2 x Ft>! — t>! x Fy 2 ). 

On the other hand, 

T 2 EF(vi,v 2 ) = vec _1 |EF(vi x v 2 )}, 
and hence (14.121) follows from the algebraic identity 

EF(vi x v 2 ) = v 2 x Fui -ui x Fw 2 , 

which holds for any F G M. 

We may further identify the four spaces of fields in (13. 3p with the corresponding spaces of 
forms in (14.111) in a natural way: 

• (u,p) G C°°(V x K) ~ (vec- 1 u,vecp) G A°(K x V). 

• F G C°°(M) ~ lo G A X (K) given by w(v) = vec~\Fv). 

• F G C°°(M) ~ /i G A 2 (V) given by ^(vi, v 2 ) = x u 2 ). 

• (w,p) G C°°(V x K) ~ G A 3 (K x V) given by u(vi,v 2 ,v 3 ) = p{v x x v 2 ■ v 3 ), 

K v h v 2, V3) = u(vi X V 2 ■ V 3 ). 

Under these identifications the we find that 

• d : A°(K) -> A : (K) corresponds to the row-wise gradient C°°(V) -> C°°(M). 

• S : A°(V) -> A X (K) corresponds to the inclusion of C°°(K) -> C°°(M). 

• d x o 5f 1 o d x : A X (K) -> A 2 (V) corresponds to J = curl S" 1 curl : C°°(M) -> C°°(M). 

• d 2 : A 2 (V) — > A 3 (V) corresponds to the row- wise divergence C°°(M) — > C°°(V). 

• 5 2 : A 2 (V) -> A 3 (K) corresponds to the operator -2skw : L7°°(M) -> L7°°(K). 
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Thus, modulo these identifications and the (unimportant) constant factor in the last identi- 
fication, (13. 3p and (14.111) are identical. Hence we have established the following result. 

Theorem 4.1. When the de Rham complex (12.21) is exact, (i.e., the domain is contractible) , 
then so is the elasticity complex (13.31) . 

To end this section, we return to the operator S : A fc (V) — > A fc+1 (K) defined by S = 
dK — Kd. Let K' : A fc (IK) — > A fc (V) be the adjoint of K (with respect to the Euclidean inner 
product on V and the Frobenius inner product on K), which is given by (K'uj) x (vi, . . . , v k ) = 
-2u x (v u v k )x. Define S' : A k (K) -> A fe+1 (V) by S' = dK' - K'd. Recall that the wedge 
product A : A fc x A' — > A k+l is given by 



(uAfi)(v 1} . . .,v k+t ) = (signer) a; 



&k+l ' 



' v °k+l )1 



E A k , jj, E A l ,Vi E V, 



where the sum is over the set of all permutations of {1, . . . , k + 1}, for which o\ < Oi < ■ ■ ■ < 
o~j. and crk+i < °~j+2 < ■ ■ • < Cfc+z- This extends as well to differential forms with values in an 
inner product space, using the inner product to multiply the terms inside the summation. 
Using the Leibniz rule (j4.6j) . we have 

(4.13) (Su) A/i = (-l) k uAS'fi, u E A fc (V), /i E A l (K). 

We thus have 



dKuj Afi= (-l) k+1 Kuj A dfx + diKuj A \i) 



-1 



,fe+i 



AK'dfi + d(u AK'/i), 



and 



Kdu A n = duj A K'n = (-l) k+1 u A dK'n + d{u A K'fj), 

Subtracting these two expressions gives (14.131) . 

For later reference, we note that, analogously to (14.71) . we have 

fc+i 



(4.14) (S'u)(v 1} ...,v k+1 ) = -2j2(-t) 



3+1 



v k +i)vj, u E A k (Q; 



5. The discrete construction 

In this section we derive a discrete version of the elasticity sequence by adapting the 
construction of the previous section. To carry out the construction, we will use two dis- 
cretizations of the de Rham sequence. For k = 0, 1,2,3, let A^ denote a finite-dimensional 
space of HA k for which dA k C A^ +1 , and for which there exist projections Uh 
which make the following diagram commute: 



Il£ : A k 



(5.1) 



A 



A 1 



A 2 



A 2 



A 3 



A 3 











This is simply the diagram (12.31) written in the language of differential forms. We do not 
make a specific choice of the discretization yet, but, as recalled in § [2j there exist many 
such discrete de Rham complexes based on piecewise polynomials. In fact, as explained in 
[7j, for each polynomial degree r > we may choose A| to be the space of all piecewise 
polynomial 3-forms with respect to some simplicial decomposition of Q, and construct four 
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such diagrams. We make the assumption that Vi{ST) C A°, which is true in all the cases 
mentioned. 

Let be a second set of finite dimensional spaces with corresponding projection operators 
lih enjoying the same properties, giving us a second discretization of the de Rham sequence. 
Supposing a compatibility condition between these two discretizations, which we describe 
below, we shall construct a discrete elasticity complex. 

We start with the complex 

(5.2) IxV^ A°(K) x A° h (V) 4- • • • 4 A£(K) x A 3 h (V) -> 

where A|(K) denotes the K- valued analogue of A^ and similarly for A*(V). For brevity, we 
henceforth write A^(W) for A*(K) x Aj°(V). As a discrete analogue of the operator K, we 
define Kh : A*(V) — > A*(K) by Kh = HhK where Uh is the interpolation operator onto 
A*(K). 

Next define S h = S kA : A|(V) -> A^ +1 (K) by 5 fc = cZA^ - A^d, for jfe = 0, 1,2. Observe 
that the discrete version of (14. 8p . 

(5.3) dS h = -S h d, 

follows exactly as in the continuous case. From the commutative diagram (15.11) . we see that 

S h = dU h K - U h Kd = Tl h (dK - Kd) = Tl h S. 

Continuing to mimic the continuous case, we define the automorphism $^ on A^(W) by 

® h (io,n) = {uj + K h fi,fj), 

and the operator A h : Af (W) -> Aj; +1 (W) by A h = fchdfcj" 1 , which leads to 

A h {u, (jl) = (dw - S h fJL, dfi). 

Thanks to the assumption that V\ C A°, we have $^(W) = $(W). Hence, inserting the 
isomorphisms $^ into (15. 2p . we obtain 

(5.4) $(W)^A°(W) ^ A£(W) ^ A|(W) ^ A£(W) -> 0. 
In analogy to the continuous case, we define 

rj = { ( w , /i) G Ai(W) | du; = S 1>h » }, F^ = { (u>, /i) G AJ(W) | w = }. 

As in the continuous case, we can identify T\ with A^(V), but, unlike in the continuous case, 
we cannot identify with A^(K), since we do not require that S^h be invertible (it is not 
in the applications). Hence, in order to derive the analogue of the diagram (I4.10p we require 
a surjectivity assumption: 

(5.5) The operator Si t h maps A ft (V) onto A^(K). 

Under this assumption, the operator Sh = S^h has a right inverse mapping A)[(K) into 
A\(Y). This allows us to define discrete counterparts of the projection operators n 1 and 7r 2 
by 

nttw, n) = (u,fi- S ] h S h n + Sldw), t^(u, //) = (0, // + dS\u), 
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and obtain the discrete analogue of ( 14. 10[) : 



$(W)^A°(W) ^ A\(W) ^ A|(W) ^ A£(W) 



-4,, 



(5.6) 



■«2 



id 



$(W)^A°(W) 



~4„ 



F 1 
1 /( 



F 2 



-4 h 



A|(W) 







It is straightforward to check that this diagram commutes. For example, if (u,fj) G A°(W), 
then 

n\A h {u, n) = 7r^(cL; - S h ji, dfi) = (du - S h fj,, d/j - SlS h dfj, + S\d[du - S h fi\) 
= (du - S h fi,dfi - Sl\S h d[i + dS h /j]) = A h (u,ii), 



where the last equality follows from (15. 3p . Thus the bottom row of (15.61) is a subcomplex of 
the top row, and the vertical maps are commuting projections. In particular, when the top 
row is exact, so is the bottom. Thus we have established the following result. 

Theorem 5.1. For k = 0, . . . , 3 ; let A^ be a finite dimensional subspace of HA k for which 
dA\ C A^ +1 and for which there exist projections lih = n£ : A k — > A\ that make the diagram 
(15. ip commute. Let A^ be a second set of finite dimensional spaces with corresponding pro- 
jection operators 11^ with the same properties. If S\ t h := dU\K — Tl 2 h Kd maps A^(V) onto 
A 2 h (K), and the bottom row of (15. ip is exact for both sequences A^ and A k , then the discrete 
elasticity sequence given by the bottom row of (15. 6p is also exact. 



The exactness of the bottom row of (15.61) suggests that the following choice of finite element 
spaces will lead to a stable discretization of (12. ip : 

£ fc ~A£(V), V h ~AftV), Q h ~A 3 h (K). 

In the next section we will make specific choices for the discrete de Rham complexes, and 
then verify the stability in the following section. 

For use in the next section, we state the following result, giving a sufficient condition for 
the key requirement that Si h be surjective. 



Proposition 5.2. If the diagram 



(5.7) 



commutes, then Si^ is surjective. 



A\N) 



A£(V) 



Si 



A 2 (K) 



A|(K) 



6. A FAMILY OF DISCRETE ELASTICITY COMPLEXES 

In this section, we present a family of examples of the general discrete construction pre- 
sented in the previous section, by choosing specific discrete de Rham complexes. These 
furnish a family of discrete elasticity complexes, indexed by an integer degree r > 0. In the 
next section we use these complexes to derive finite elements for elasticity. In the lowest 
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order case, the method will require only piecewise linear functions to approximate stresses 
and piecewise constants to approximate the displacements and multipliers. 

We begin by recalling the two principal families of piecewise polynomial spaces of differ- 
ential forms, following the presentation in [7J. We henceforth assume that the domain Q is 
a contractible polyhedron. Let %i be a triangulation of fl Let Th be a triangulation of Q by 
tetrahedra, and set 

V r A k {%) = {uoe HA k (Q) | u\ T G V r A k (T) VTeT h }, 
V+A k (T h ) = {ue HA k (Q) | u\ T G V+A k {T) WTe%}. 

Here V+A k (T) := V r A k (T) + nV r A k+1 (T) where k : A fc+1 (T) -> A k (T) is the Koszul differ- 
ential defined by 

{kuj) x {v 1 , ■ ■ ■ , v k ) = uj x (x, v 1 ,--- ,v k ). 

The spaces VfA°(Th) = V r +iA° (Th) correspond to the usual degree r + 1 Lagrange piecewise 
polynomial subspaces of H 1 , and the spaces VfA 3 (Th) = V r A 3 (Th) correspond to the usual 
degree r subspace of discontinuous piecewise polynomials in L 2 (Q). For k = 1 and 2, 
the spaces V^A k {Th) correspond to the discretizations of if (curl) and H(div), respectively, 
presented by Nedelec in [25], and the spaces V r A k (7h) are the ones presented by Nedelec in 
[25] . An element uj G V r A k (T h ) is uniquely determined by the following quantities: 

(6.1) J^uAC, (6^ r VwA"(/), feA d (T h ), k<d<3. 

Here A d (T h ) is the set of vertices, edges, faces, or tetrahedra in the mesh T h , for d — 0,1, 2, 3, 
respectively, and for r < 0, we interpret VfA k (T) = V r A k (T) = 0. Note that for u G A k , 
u naturally restricts on the face / to an element of A k (f). Therefore, for C G A d ~ k (f), the 
wedge product uj A ( belongs to A d (f) and hence the integral of uj A £ on the ci-dimensional 
face / of T is a well-defined and natural quantity. Using the quantities in (16.11) . we obtain a 
projection operator from A fe to V r A k (Th). 

Similarly, an element uj G V^A k (T h ) is uniquely determined by 

(6.2) 1 UAC Ce^A"(/), / G A d (T h ), k<d<3, 

and so these quantities determine a projection. 

If X is a vector space, we use the notation P r A fc (7/ 1 ;X) and VfA k (Th; X) to denote the 
corresponding spaces of piecewise polynomial differential forms with values in X. Further- 
more, if X is an inner product space, the corresponding degrees of freedom are given by (16.11) 
and (16.21) . but where the test spaces are replaced by the corresponding X valued spaces. 

To carry out the construction described in the previous section we need to choose the two 
sets of spaces A k and A| for k = 0,1, 2, 3. We fix r > and set A k = V+A k (T h ), k = 0,1, 2, 3, 
and A° h = V r+2 A (T h ), A\ = V+ +1 A l {T h ), A\ = V r+l A 2 {T h ), and A\ = P r A\. As explained in 
[7J, both these choices give a discrete de Rham sequence with commuting projections, i.e., a 
diagram like (15. ip makes sense and is commutative. 

We establish the key surjectivity assumption for our choice of spaces, by verifying the 
commutativity of (15.71) . 
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Lemma 6.1. Let A*(V) = V+ +1 A\T h ;Y) and A 2 h {K) = V+A 2 (T h ;K) with projections TL h , 
U\ defined via the corresponding vector valued moments of the form (16. ip and (I6.2p . If 
S^h = then 

(6.3) S 1>h Il 1 h = n 2 h s u 
and so Sx,h is surjective. 

Proof. We must show that (U 2 h Si - S^U^a = for a E A X (V). Defining u — (I — ft^cr, 
the required condition becomes U^Siou = 0. Since U\uj = 0, we have 

(6.4) jfwAC = 0, Ce^A^^Y), / G A d (T h ), 2 < d < 3, 

(in fact (16. 4p holds for d = 1 as well, but this is not used here). We must show that (16.41) 
implies 

(6.5) Js 1 uA fJ , = 0, ^P r+2 _/- 2 (/;K), / G A d (T h ), 2<d<3. 

From (Q31) . we have SiU) A /i = -w A £ where C = S^A* e V r+ 2-dA d - l {f; V) for /i G 
P r+2 _dA a! - 2 (/;K), as is evident from (|P3)) . Hence O follows from (JS3D- □ 

7. Stable mixed finite elements for elasticity 

Based on the discrete elasticity complexes just constructed, we obtain mixed finite element 
spaces for the formulation (12.11) of the elasticity equations by choosing E^, Vh, and Qh as 
the spaces of matrix and vector fields corresponding to appropriate spaces of forms in the 
K- and V-valued de Rham sequences used in the construction. Specifically, these are 

(7.1) Z h ~V r+1 A 2 (T h ;Y), V h ~p r A*(T h ;Y), Q h ~ V r A 3 (T h ; K). 

In other terminology, E^ may be thought of as the product of three copies of the Nedelec 
H(div) space of the second kind of degree r + 1, and Vh and Qh are spaces of all piecewise 
polynomials of degree at most r with values in K and V, respectively, with no imposed 
interelement continuity. In this section, we establish stability and convergence for this finite 
element method. 

The stability of the method requires the two conditions (Al) and (A2) stated in § 2. 
The first condition is obvious since, by construction, divE ft C Vh, i.e., dV r+ iA 2 (T h ;Y) C 
V r A 3 (Th,Y). The condition (A2) is more subtle. We will prove a stronger version, namely 

(A2') for all nonzero (v, q) E Vh X Qh, there exists nonzero r E E^ with divr = v, 
2IlQ h skw t = q and 

l|r||div<c(|H| + y|), 

where Hq h is the L? projection into Qh and c is a constant. 

Recalling that = x T , r+ iA 2 {T h ;\) and Ah(0,a) = (— S^cr, da), and that the operator 
S2 corresponds on the matrix level to —2 skw, we restate condition (A2') in the language of 
differential forms. 
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Theorem 7.1. Given (u,fi) G V r A 3 (T h ;K) x V r A 3 (T h ; Y), there exists a G V r+ iA 2 (T h ;Y) 
such that Ah(0, a) = (u, fi) and 

(7-2) H HA <c(||u,|| + |M|), 

where the constant c is independent ofu,fx and h. 

Before proceeding to the proof, we need to establish some bounds on projection operators. 
We do this for the corresponding scalar-valued spaces. The extensions to vector-valued 
spaces are straightforward. First we claim that 

(7.3) \\U 2 h ri\\ < c|M|i V^Gi^A 2 , ||IT^;|| < c||w|| G i^A 3 . 

Here the constant may depend on the shape regularity of the mesh, but not on the meshsize. 
The second bound is obvious (with c = 1), since 11^ is just the L 2 projection. The first 
bound follows by a standard scaling argument. Namely, let T denote the reference simplex. 
For any (3 G V r+ iA 2 (T), we have 

(7.4) Piio,f <^(EEi /^ a ai + Ei /^ a cd> 

/ A f c T 

where / ranges over the faces of T, fi over a basis for V^~A°(f), and ( over a basis for 
Vr-iA l (T). This is true because the integrals on the right hand side of (I7.4p form a set of 
degrees of freedom for j3 G V r+ iA 2 (T) (see f l6.ll) ). and so we may use the equivalence of all 
norms on this finite dimensional space. We apply this result with {3 = tl^fj, where tl\ is the 
projection defined to preserve the integrals on the right hand side of (17. 4p . It follows that 

lin^llo.f <c(EEl /^ A ^ + Z)l /^ A CI) < 4v\\i,f 



Ti 



where we have used a standard trace inequality in the last step. Next, if T is an arbitrary 
simplex and 77 G if x A 2 (T), we map the reference simplex T onto T by an affine map x 1— > 
Bx + 6, and define 17 G H 1 A 2 (f) by 

mi^h = Vx(Bv\, BV2), 
for any x = Bx + b G T and any vectors Vi,V2- It is easy to check that fl^r] = fl^f), and that 
llfifc^llo.r < c||n 2 ?)|| 0> f < c\\rj\\ lt f < c(||?7||o,t + h\rj\ 1:T ) < c||t7|| 1)T . 

Squaring and adding this over all the simplices in the mesh 7/j gives the first bound in (17. 3p . 

We also need a bound on the projection of a form in H x A l into A^ = 'P r f fl A 1 (7/ l ). However, 
the projection operator 11^ is not bounded on H 1 , because its definition involves integrals 
over edges. A similar problem has arisen before (see, e.g., [10J ) , and we use the same remedy. 
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Namely we start by defining an operator : i/ 1 A 1 — > V^ + iA 1 (%) by the conditions 

(7.5) / til h u A C = / u A C, C e V r ^A 2 (T), T G T h , 
Jt Jt 

(7.6) Jti$ h uA{ = JwA(, (EV r A\f), feA 2 (T h ), 

(7.7) y*nJ h o;AC = J Cen+iA°(e), eeA^). 

Note that, in contrast to 11^, in the definition of iL 1 ^, we have set the troublesome edge 
degrees of freedom to zero. Let : i^A^T) -> P r + ^(T) be defined analogously on the 
reference element. 

Now for p e ^A 1 ^), dUlp E V r +iA 2 {T ), so 

||dn^|| 0>f <c(EEl /^oPAA| + X)| [dfllpACD, 

f fi f c T 

where again / ranges over the faces of T, /i over a basis of P+A°(/), and £ over a basis of 
V^_ 1 A 1 {f). But 



^ dfloP A p — J LTqP A dp = J p A dp, 



where we have used Stokes theorem and the fact that the vanishing of the edge quantities 
in the definition of LTq to obtain the first equality, and the face degrees of freedom entering 
the definition of LTq to obtain the second. Similarly, 

[dfllp A( = ftllpAdl+f fl l p AC = [jAd(+ [ pAC= [dp AC 
Jt Jt J&t Jt J&t Jt 

It follows that 

||«fll5p|| ,T < 4f>\\i,t , pe^A'if). 

When we scale this result to an arbitrary simplex and add over the mesh, we obtain 

WdU^pW ^c^HpH + IIpIK), peH'A 1 ^). 

To remove the problematic hr x in the last estimate, we introduce the Clement interpolant 
Rh mapping if 1 A 1 into continuous piecewise linear 1-forms (still following [10]). Then 

||p - R hP \\ + h\\ P - R hP \\ x < ch\\ P \\ u P e H l A x . 

Defining TL\ : H 1 A 1 -> V+ +1 A\ by 

(7-8) til = TLl h (I-R h ) + R h , 

we obtain 

\\du\pW < Hrffi^cx-^pii + ii^^ii < cC^ikz-^)^!!^!!^-^)^!!!^!!^^!!) < cUpIK- 

Thus we have shown that 

(7.9) \\<mlp\\ < cWpWu peH x K\ 
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Having modified 11^ to obtain the bounded operator LL^, we now verify that the key 
property (16.31) in Lemma [6.11 carries over to 

(7.io) s 1A ni = n^i, 

where we now use the vector-valued forms of the projection operators. It follows easily from 
(EU), CHD, and flUD that (EHJ) holds with uj = (I - Ti\)a, so that the proof of fl7A0|) is the 
same as for (16.31) . 

We can now give the proof of Theorem 17.11 

Proof of Theorem 7.1, Given p G 7 :> r A 3 (7/ l ;V) there exists n G H 1 A 2 (Y) such that dr] = 
p with the bound ||r/||i < c\\/jl\\ (since d maps H l A 2 onto L 2 A 3 ). Similarly, given uj G 
•p r A 3 (7^;K) there exists r G H l A 2 (K) with dr = uj + S 2>h fllr] with the bound ||r||i < 
c\\co + S^htl^vW- Let p = S^t (recall that Si is an isomorphism) and set 

a = dUlp + tllrj. 

We will now show that .4.^(0, a) = (uj, p). From the definition of a, we have 

S 2 ,hO- = -S 2 ,hdfllp - S 2 ,hfl 2 h V- 
Then, using (15. 3p . (17.101) . and the commutativity dllh = Uhd, we see 

S 2 ,hdtLlp = -dS ljh TLlp = -dUlSxp 

= -dU 2 r = -U 3 h dr = -U 3 h (u + S 2 , h fl 2 hV ) = -uj- S 2 , h fl 2 hV . 

Combining, we get —S 2 ^a = to as desired. Further, from the commutativity dllh = Tlhd and 
the definition of n, we get 

da = dXl 2 h r) = fl 3 h drj = flip = p, 

and so we have established that Ah(0,a) = (p,uj). 

It remains to prove the bound (17.21) . Using (17. 3p . we have 

H^n^n = lin^^n < c\\s 2 fi 2 hV \\ < c ||n^|| < chh < c\\p\\. 

Thus lipid < c||r||i < c(|H| + ||/x||). Using (EHJ) , we then get \\dn.\p\\ < c||p||i < c(|M| + |M|), 
and, using (17.31) . that ||n^|| < c IM|i < c l I A*II - Therefore \\a\\ < c(\\uj\\ + \\p\\), while 
1 1 da 1 1 = \\p ||, and thus we have the desired bound (17.21) . □ 

We have thus verified the stability conditions (Al) and (A2), and so may apply the stan- 
dard theory of mixed methods (cf. [12], [13], p2], [20]) and standard results about approxi- 
mation by finite element spaces to obtain convergence and error estimates. 

Theorem 7.2. Suppose (a,u,p) is the solution of the elasticity system (11.51) and (ah,Uh,Ph) 
is the solution of discrete system (12.11) . where the finite element spaces S/j, Vh, and Qh are 
given by (17.11) for some integer r > 0. Then there is a constant C, independent of h, such 
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that 

\W - cr h \\ div + \\u-u h \\ + \\p-Ph\\ < C inf (||cr - r h || div + + ||p-g&||), 

T h EY, h ,v h €V h ,q h eQ h 

h -*h\\ + lb - pall + Ik - n>|| < c(\\a - nrV|| + |b - rr>||), 
lb - « fc || <c(\\a- nrV|| + lb - Kp\\ + \\u- fq u \\), 

\\d(a - a h )\\ = \\da-IL n da\\. 
If u and a are sufficiently smooth, then 
\\a — 0^11 + |b - u h \\ + |b - Ph\\ < C , /i r+1 |b|| r+2 , || div(a - a h )\\ < Ch r+1 \\ div cr|| r+1 . 

Remark. Note that the errors ||a — o~h\\ an d \\iih — n^u|| depend on the approximation of 
both a and p. For the choices made here, the approximation of p is one order less than the 
approximation of a, and thus we do not obtain improved estimates, as one does in the case 
of the approximation of Poisson's equation, where the extra variable p does not enter. 

8. A SIMPLIFIED ELEMENT 

Recall that the lowest order element in the stable family described above, for a discretiza- 
tion based on (11.51) . is of the form 

E h ~ V x t?(J h - V), V h ~ V A 3 (T h ; V), Q h ~ V Q A 3 (T h ;K). 

The purpose of this section is to present a stable element which is slightly simpler than this 
one. More precisely, the spaces Vh and Qh are unchanged, but will be simplified from full 
linears to matrix fields whose tangential-normal components on each two dimensional face 
of a tetrahedron are only a reduced space of linears. 

We will still adopt the notation of differential forms. By examining the proof of Theo- 
rem 17.11 we realize that we do not use the complete sequence (15. 2p for the given spaces. We 
only use the sequences 

(8 1} V+A 2 (T h] K) 4. V A 3 (T h] K) - 0, 

V+A\T h] N) 4 V x A 2 (T h] N) 4 P A 3 (T,;V) - 0. 

The purpose here is to show that it is possible to choose subspaces of some of the spaces in 
(18. ip such that the desired properties still hold. More precisely, compared to (18. ip . the spaces 
P ] f A 1 (7/ l ;V) and PiA 2 (7^;V) are simplified, while the three others remain unchanged. If 
we denote by Vt-A l {T h ;N) and V 1 -A 2 (T h ;Y) the simplifications of the spaces V^A 1 (T h ;Y) 
and ViA 2 (Th,Y), respectively, then the properties we need are that: 

(8-2) V^A\T h ;Y) 4 7VA 2 (T,;V) 4 V A 3 (T h] Y) - 

is a complex and that the surjectivity assumption (15.51) holds, i.e., Sh = 5*1, h maps the space 
Vf_A\T h ;Y) onto V^A 2 {T h ;K). Note that if P +A 2 (T h ; V) C P l5 _A 2 (T h ;V), then d maps 
Vi,-A 2 (T h ;Y) onto V A 3 {T h ; V). 

We first show how 7 ? 1 f _A 1 (7/ l ; V) can be constructed as a subspace of Vi A 1 {T h ;Y) . Since 
the construction is done locally on each tetrahedron, we will show how to construct a space 
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"P^A^T; V) as a subspace of P ] f A 1 (T; V). We begin by recalling that the face degrees of 
freedom of A X (T; V) have the form 

[uAfj,, fj,eV A l {f,V). 
Jf 

We then observe that this six dimensional space can be decomposed into 

V^if- V) = V A l {f- T f ) + V A\f; N f ), 

i.e. into forms with values in the tangent space to /, Tf or the normal space Nf. This is a 
4 + 2 dimensional decomposition. Furthermore, 

V A\f; Tf) = V Al ym (f; T f ) + V^{J\ 

where \i G 7 ? A 1 (/; Tf) is in "PoAg (/; Tf) if and only if fi(s) ■ t = fi(t) ■ s for orthonormal 
tangent vectors s and t. Finally, we obtain a 3 + 3 dimensional decomposition 

V A\f; V) = V Al ym (f; T f ) + VoA^f; V), 

where 

P A s 1 kw (/; V) = V Al kw (f; T f ) + V A\f; Nf). 
In more explicit terms, if \x G VqA 1 (F] V) has the form 

H(q) = (ait + a 2 s + a^n)q ■ t + (a^t + a$s + a§n)q ■ s, 

where t and s are orthonormal tangent vectors on /, n is the unit normal to /, and q 
is a tangent vector on /, then we can write u = /ii + /i 2 , with \i\ G T^Agy^/; V) and 
/^GPoA^CAV), where 

/ a 2 + a 4 \ /a 2 + a 4 \ 
yUi(g) = I a x t H s 1 g ■ t + I 1 + a 5 sjq-s, 

^2{q) = I ^ s + a 3n) q-t + i 1 + a 6 njq-s. 

The reason for this particular decomposition of the degrees of freedom is that if we examine 
the proof of Lemma 16.11 where equation (16.31) is established, we see that the only degrees of 
freedom that are used for an element u> G A 1 (T; V) are the subset of the face degrees of 
freedom given by: 

JuA(S' u), veV A°(f;K). 

However, for v G 7 : oA (/;lK), fj, = S' u is given by «(g) = vq. Since the general element 
v G 7 ? oA°(lK) can be written in the form bi(ts T — st T ) + b2(nt T — tn T ) + b 3 (ns T — sn T ), 
vq = (—bis + b 2 n)q ■ t + {pit + b 3 n)q ■ s for q a tangent vector, and thus u G 'PoAs kw (/; V). 
Hence, we have split the degrees of freedom into three on each face that we need to retain 
for the proof of Lemma 16.11 and three on each face that are not needed. The reduced space 
Vi_A l (T; V) that we now construct has two properties. The first is that it still contains the 
space PiA : (T; V) and the second is that the unused face degrees of freedom are eliminated 
(by setting them equal to zero). We can achieve these conditions by first writing an element 
uj G Vi A l {T; V) as u> = Uf l u>+(I— Hh)u, where Uh denotes the usual projection operator into 
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'PiA 1 (T; V) defined by the edge degrees of freedom. Then the elements in (J— Xlh)Vi A X (T; V) 
will satisfy 

uAfi = 0, fj, G PiA°(e; V), e G Ai(T), 

i.e., their traces on the edges will be zero. Thus, they are completely defined by the face 
degrees of freedom 

JuAfi, neVoA^fiV), feA 2 (T). 

Since this is the case, we henceforth denote (J — Hh)Vi A : (T; V) by P 1 f yA 1 (T; V). 
We then define our reduced space 

Pf_A x (T;V) = T^i A 1 (T; V) +P+ /j _A 1 (T; V), 
where P^_A 1 (T; V) denotes the set of forms to G P-j^A 1 (T; V) satisfying 

[uA f i = 1 /ieP A^(/;¥), 

i.e., we have set the unused degrees of freedom to be zero. 
Then 

V+_Al(T h ;Y) = {ue V+A\T h ;Y) : u\ T G P+.A 1 ^; V), VT G TJ. 
The degrees of freedom for this space are then given by 

(8.3) JuAfJi, \i G "PiA (e;V),e G Ai(T), jfwAfx, /x 6 V A x 8kw {f\ V), / G A 2 (T). 

It is clear from this definition that the space P ] f _A 1 (T; V) will have 48 degrees of freedom 
(36 edge degrees of freedom and 12 face degrees of freedom). The unisolvency of this space 
follows immediately from the unisolvency of the spaces ViA l (T; V) and _A 1 (T; V). 

The motivation for this choice of the space Vi_A\(Th,Y) is that it easily leads to a 
definition of the space Vi-A 2 (T h ;Y) that satisfies the property that (18. 2p is a complex. We 
begin by defining 

?VA 2 (T; V) = P +A 2 (T; V) + dV^A^T; V). 

It is easy to see that this space will have 24 face degrees of freedom. Note this is a reduction 
of the space ViA 2 (T; V), since 

ViA 2 {T; V) = P + A 2 (T; V) + dVl f A\T; V). 

We then define 

7VA 2 (T ft ;V) = {u G V 1 A 2 (T h ; V) : u\ T G ?VA 2 (T; V), VT G TJ. 

It is clear that V^A 2 (T h ;V) C P lt -A 2 (T h ;V). The fact that the complex is exact now 
follows directly from the fact that the complex 

(8.4) PiA^TjV) 4 P+A 2 (T;V) 4 P A 3 (T;V) -> 
is exact and the definition 

Pl-A^TjV) =P + A 1 (T;V) + P+ /j _A 1 (T;V). 
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We will define appropriate degrees of freedom for the space Vi-A 2 (T; V) by using a subset 
of the 36 degrees of freedom for V\A 2 (T; V), i.e., of fj uj A fi, [i G ViA°(f; V). In particular, 
we take as degrees of freedom for Vi-A 2 (T; V), 

JuAn, neVi^A^y), V/eA 2 (T), 

where V liS k w A (f;V) denotes the set of fi G ViA°(f; V) that satisfy d/i G V A\ kw {f; V). It 
is easy to check that such will have the form 

(8.5) fj, = fiQ + • £)n + 0:2(2; • s)n + 0:3 [(a; • t)s — (x ■ s)t], 

where fi G V A°(f;Y). 

Since Vi tS k w A°(f;Y) is a 6-dimensional space on each face, the above quantities specify 
24 degrees of freedom for the space Pi i _A 2 (T;V). To see that these are a unisolvent set 
of degrees of freedom for P li _A 2 (T; V), we let uj = u + dui, where u G VqA 2 (T; V) and 
LJi G Vij-A 1 ^; V) and set all degrees of freedom equal to zero. Then for \i G P A°(/; V), 
since 

J {ujo + dui) A /i = J u) A fj,, 

we see that c^o = by the unisolvency of the standard degrees of freedom for VqA 2 (T; V). 
In addition, for all fi G Vi tS k w A°(f; V) and uq = 0, 



w A [A — / da;i A (J, — I uj\ A <i/x. 
/ J f J f 

Since d/i G VoAl kw (f;Y), oj\ = by the unisolvency of the degrees of freedom of the space 
Vx, f ,-A l (T;V). 

Using an argument completely parallel to that used previously, it is straightforward to 
show that the simplified spaces also satisfy assumption (15.51) . i.e., that Sh is onto. 

To translate the degrees of freedom of the space Vi t -A 2 (T;Y) to more standard finite 
element degrees of freedom, we use the identification of an element u G A 2 (T;V) with a 
matrix F given by u(vi,V2) = F{y\ x t> 2 ). Then to(t,s) = Fn and J^w A/i = J fi T Fndf. 
Since /i G Vi tS k w A°(f] V) and hence is of the form (|8.5p . we get on each face the six degrees 
of freedom 

J^Fndf, J(x -t)n T Fndf, J (x ■ s)n T Fndf, J \(x ■ t)s T - {x ■ s)t T ]Fndf. 

Finally, we note that the analogue of Theorem (17. 2p holds with r = for the simplified 
spaces. 
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